function newquadtest()
clf
pitch = 1.0;
number_points = 27;
number_angles = 2*floor(log3(number_points))+1;

% generate unique sines
sines = zeros(number_angles, 1);
N = ceil(number_angles/2);
for i = 1:N
    if i == 1
        num = number_points;
    else
        num = 3 ^ (N-i);
    end
    den = sqrt( num^2 + 1 );
    sines(i)=num/den;
end
for i = 1:N-1
    if i == 1
        num = N-1;
    else
        num = 3 ^ (N-i);
    end
    sines(i+N) = 1/sqrt( 3^(2*i) + 1 );
end
sines
phi = asin(sines)*180/pi
w = weights(phi), sum(w)

p = (asin(sines)*180/pi-45)/45;
delta(1) = 1-p(1);
for i = 2:length(p)
    delta(i) = p(i-1)-p(i)-delta(i-1);
end
delta'
sum(delta)
% number_angles
% maxx=number_angles-1;
% xx = 0:maxx;
% phicheb = ((45 + 45*cos(xx * pi/(maxx))) - 45)*84.8056/90 + 45
% generate tracks originate
coses = sqrt(1-sines.^2);
s = linspace(0,1,100);
c = sqrt(1-s.^2);
% ss = [     9.876883405951378e-01
%      8.910065241883678e-01
%      7.071067811865475e-01
%      4.539904997395467e-01
%      1.564344650402309e-01];
% cc = sqrt(1-ss.^2);
xx = linspace(0,1,10);
plot(c,s,'k--',coses,sines,'b.','MarkerSize',16,'LineWidth',2)
hold on
for i = 1:10
   plot( [xx(i) xx(i)], [0 1], 'k--')
   plot( [0 1],[xx(i) xx(i)], 'k--')
end

axis square




end

function y=log3(x)
    y = log(x)/log(3); 
end

function w=weights(phi)
disp('weights')
phi = phi/90
n = length(phi);
u = zeros(n,1);
A = zeros(n);
for i = 1:n
    % compute analytic u_i = int(0,1,x^i)
    u(i) = 1/i;
    % row of A
    A(i,:) = phi .^ (i-1);
end
w = A\u;
end

function plotterCart(track, MeshX, MeshY, Index)
% function plotter(track, Pitch, Radii, Index)
% 
% Input:
%       track       - tracking data
%       MeshX       - Spatial divisions along x axis
%       MeshY       - Spatial division along y axis
%       Index       - angular index ( from 1 to NumAzimuth )
% Output:
%       a plot
%

    Pitch = MeshX(end)-MeshX(1);
    
    % Plotting square
    Xa = [ 0.0; Pitch; Pitch; 0.0; 0.0 ];
    Ya = [ 0.0; 0.0; Pitch; Pitch; 0.0 ];
    plot (Xa ,Ya ,'k','LineWidth',4);
    hold on;
    
    % Plotting meshes
    k = 0;
    for i = 1:length(MeshX)-1
        for j = 1:length(MeshY)-1
            k = k+1;
            x1 = MeshX(i); x2 = MeshX(i+1);
            y1 = MeshY(j); y2 = MeshY(j+1);
            X = [x2 x2 x1 x1]; Y = [y2 y1 y1 y2];
            fill(X,Y,col(k));        
        end
    end
    
    % Plotting tracks
    I = squeeze( track.Enter(Index,:,:) );
    F = squeeze( track.Exit(Index,:,:) );
    for k = 1:track.NumSpace
        plot ([I(k ,1) ;F(k ,1) ],[I(k ,2) ;F(k ,2) ],'g','LineWidth',2);
    end
    
    axis ([0 Pitch 0 Pitch])
    axis square ;
    grid on;
    hold off;
    
    xlabel(' x [cm] ')
    ylabel(' y [cm] ')
    title([' Pin Cell Tracking -- Angle Index ', num2str(Index)])
end

function color = col(g)
    % this function is a hard-coded color map for the different
    % flux groups.  I've accounted for up to 8 groups.
    % set(ur,'Color',[1 0.7 0.2],'LineWidth',2);
    switch g
        case 1
            color = [0.0 0.0 1.0]; % blue
        case 2
            color = [0.0 0.8 0.2]; % nice green
        case 3
            color = [1.0 0.0 0.0]; % red
        case 4
            color = [0.4 0.0 0.6]; % purple
        case 5
            color = [0.9 0.4 0.0]; % orange
        case 6
            color = [0.5 0.2 0.0]; % brown
        case 7
            color = [0.0 0.8 0.6]; % turquoise
        case 8
            color = [0.7 0.6 0.0]; % gold
        otherwise
            color = [rand rand rand]; % random
    end
end
